Python with PostgreSQL & PostGIS¶
Note: Please always run the complete Jupyter Notebook from the beginning, as object names such as 'sql' and 'gdf' are reused in the code cells.
Libraries and Settings¶
The 3 types of swiss coordinate systems are
We adapt the code in this part: ST_transform(p.way, 4326) as geom
We have either 21781, 20256 or 4326
This where I can find the information regarding shops code
https://wiki.openstreetmap.org/wiki/Key:shop
TO get to coordinates https://tools.retorte.ch/map/
# Libraries
import os
import folium
import pandas as pd
import geopandas as gpd
from sqlalchemy import create_engine, text
# Ignore warnings
import warnings
warnings.filterwarnings("ignore")
print(os.getcwd())
/workspaces/python_postgresql_postgis
Create database connection¶
# Set up database connection
user = "pgadmin"
password = "geheim"
host = "localhost"
port = "5432"
database = "osm_switzerland"
# Create Connection URL
db_connection_url = "postgresql://" + user + ":" + password +\
"@" + host + ":" + port + "/" + database
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Test database connection
with engine.connect() as connection:
result = connection.execute(text('SELECT current_database()'))
print(result.fetchone())
# Dispose the engine
engine.dispose()
('osm_switzerland',)
List tables in database¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Open a connection
with engine.connect() as connection:
# Execute the query
result = connection.execute(text("""SELECT table_name
FROM information_schema.tables
WHERE table_schema = 'public';"""))
# Fetch and print the results
for row in result:
print(row[0])
# Dispose the engine
engine.dispose()
geography_columns geometry_columns spatial_ref_sys planet_osm_roads planet_osm_line planet_osm_polygon planet_osm_point
Show columns and data types of selected table¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Specify your table name
table_name = 'planet_osm_polygon'
# Query to get column information
query = f"""SELECT column_name, data_type
FROM information_schema.columns
WHERE table_name = '{table_name}'"""
# Execute the query and read the result into a DataFrame
df = pd.read_sql(query, engine)
# Dispose the engine
engine.dispose()
# Print the DataFrame
df
| column_name | data_type | |
|---|---|---|
| 0 | osm_id | bigint |
| 1 | z_order | integer |
| 2 | way_area | real |
| 3 | way | USER-DEFINED |
| 4 | addr:housenumber | text |
| ... | ... | ... |
| 68 | wood | text |
| 69 | tracktype | text |
| 70 | access | text |
| 71 | addr:housename | text |
| 72 | addr:street | text |
73 rows × 2 columns
Query: Select buildings for which full address is available in defined zip code areas¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query
sql = """SELECT
p.osm_id,
p."addr:street",
p."addr:housenumber",
p."addr:city",
p."addr:postcode",
p.building,
st_transform(p.way, 4326) AS geom
FROM
public.planet_osm_polygon AS p
WHERE
p."addr:street" IS NOT NULL
AND p."addr:housenumber" IS NOT NULL
AND p."addr:city" IS NOT NULL
AND p."addr:postcode" IN ('8001', '8002')"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)
# Dispose the engine
engine.dispose()
# Print the GeoDataFrame
gdf
| osm_id | addr:street | addr:housenumber | addr:city | addr:postcode | building | geom | |
|---|---|---|---|---|---|---|---|
| 0 | 39517505 | Schanzengasse | 14 | Zürich | 8001 | yes | POLYGON ((8.54989 47.36681, 8.54996 47.36676, ... |
| 1 | 288353536 | Fortunagasse | 34 | Zürich | 8001 | yes | POLYGON ((8.53967 47.37337, 8.5397 47.37334, 8... |
| 2 | 288353533 | Fortunagasse | 30 | Zürich | 8001 | yes | POLYGON ((8.53967 47.37346, 8.53969 47.37344, ... |
| 3 | 288353542 | Rennweg | 42 | Zürich | 8001 | yes | POLYGON ((8.53931 47.37351, 8.53933 47.37347, ... |
| 4 | 288353540 | Rennweg | 38 | Zürich | 8001 | yes | POLYGON ((8.53933 47.37347, 8.53939 47.37338, ... |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 2350 | 108633058 | Brandschenkestrasse | 152 | Zürich | 8002 | industrial | POLYGON ((8.52413 47.36489, 8.52415 47.36477, ... |
| 2351 | 108633096 | Brandschenkestrasse | 152c | Zürich | 8002 | office | POLYGON ((8.52388 47.36409, 8.52391 47.3639, 8... |
| 2352 | 108633055 | Brandschenkestrasse | 152b | Zürich | 8002 | industrial | POLYGON ((8.5239 47.36435, 8.5239 47.36431, 8.... |
| 2353 | 108626600 | Brandschenkestrasse | 110a | Zürich | 8002 | office | POLYGON ((8.52458 47.3655, 8.52462 47.36534, 8... |
| 2354 | 34572240 | Brandschenkestrasse | 110 | Zürich | 8002 | commercial | POLYGON ((8.52469 47.36579, 8.52476 47.36552, ... |
2355 rows × 7 columns
Show selected features on map¶
Note the popup field in the map, which has been added to provide additional information about buildings.
Example of alternative background maps (maptiles) are:
- EsriWorldImagery
- EsriWorldTopoMap
- EsriWorldGrayCanvas
- CartoDBDarkMatter
- CartoDBPositron
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=15,
tiles='EsriWorldImagery')
# Map settings
folium.GeoJson(
gdf,
name='geojson',
weight=0.5,
fill_color='greenyellow',
fillOpacity=0.8,
popup=folium.GeoJsonPopup(fields=['addr:street',
'addr:housenumber',
'addr:city',
'addr:postcode',
'building'])
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Doing the same exercise for my city: Treytorrens¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query
sql = """SELECT
p.osm_id,
p."addr:street",
p."addr:housenumber",
p."addr:city",
p."addr:postcode",
p.building,
st_transform(p.way, 4326) AS geom
FROM
public.planet_osm_polygon AS p
WHERE
p."addr:street" IS NOT NULL
AND p."addr:housenumber" IS NOT NULL
AND p."addr:city" IS NOT NULL
AND p."addr:postcode" IN ('1530')"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)
# Dispose the engine
engine.dispose()
# Print the GeoDataFrame
gdf
| osm_id | addr:street | addr:housenumber | addr:city | addr:postcode | building | geom | |
|---|---|---|---|---|---|---|---|
| 0 | 690348871 | Aéropôle | 201 | Payerne | 1530 | industrial | POLYGON ((6.91756 46.83402, 6.91788 46.83383, ... |
| 1 | 138739158 | Grand'Rue | 16 | Payerne | 1530 | yes | POLYGON ((6.93893 46.82179, 6.93904 46.82161, ... |
| 2 | 138739168 | Grand'Rue | 20 | Payerne | 1530 | yes | POLYGON ((6.93863 46.82193, 6.93868 46.82187, ... |
| 3 | 138739172 | Grand'Rue | 22 | Payerne | 1530 | yes | POLYGON ((6.93866 46.82174, 6.93872 46.82158, ... |
| 4 | 138739160 | Rue des Blanchisseuses | 18 | Payerne | 1530 | apartments | POLYGON ((6.93818 46.82203, 6.93828 46.8219, 6... |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 1045 | 357477601 | Rue des Jumelles | 22 | Payerne | 1530 | yes | POLYGON ((6.9367 46.81495, 6.9368 46.81489, 6.... |
| 1046 | 357477604 | Rue des Jumelles | 24 | Payerne | 1530 | yes | POLYGON ((6.93651 46.81481, 6.9366 46.81474, 6... |
| 1047 | 357477606 | Rue des Jumelles | 25 | Payerne | 1530 | yes | POLYGON ((6.93677 46.81465, 6.93685 46.8146, 6... |
| 1048 | 357477610 | Rue des Jumelles | 29 | Payerne | 1530 | yes | POLYGON ((6.93669 46.81428, 6.9368 46.81422, 6... |
| 1049 | 357477608 | Rue des Jumelles | 27 | Payerne | 1530 | yes | POLYGON ((6.93634 46.81433, 6.93643 46.81428, ... |
1050 rows × 7 columns
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=17,
tiles='EsriWorldTopoMap')
# Map settings
folium.GeoJson(
gdf,
name='geojson',
weight=0.5,
fill_color='red',
fillOpacity=0.9,
popup=folium.GeoJsonPopup(fields=['addr:street',
'addr:housenumber',
'addr:city',
'addr:postcode',
'building'
])
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Query: Select coffee stores in Switzerland¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query
sql = """SELECT
h.osm_id,
h.shop,
h.name,
ST_Transform(h.way, 4326) AS geom
FROM planet_osm_point h
WHERE h.shop = 'coffee';"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)
# Dispose the engine
engine.dispose()
# Print the GeoDataFrame
gdf
| osm_id | shop | name | geom | |
|---|---|---|---|---|
| 0 | 10275946038 | coffee | Quinta Coira | POINT (9.53125 46.8485) |
| 1 | 11920528580 | coffee | Maison du Cafe | POINT (8.98798 45.87081) |
| 2 | 9780147859 | coffee | Bel-Air Coffee | POINT (6.62851 46.52281) |
| 3 | 11930216622 | coffee | Charbon Café | POINT (6.6277 46.51435) |
| 4 | 7240932670 | coffee | Nespresso | POINT (6.69463 46.50053) |
| ... | ... | ... | ... | ... |
| 113 | 11993710081 | coffee | Berger Kaffee | POINT (7.62597 46.87853) |
| 114 | 1438791023 | coffee | Nespresso | POINT (7.63094 46.75724) |
| 115 | 1517037969 | coffee | Nespresso | POINT (7.15165 46.80227) |
| 116 | 5019777080 | coffee | KAFOJ | POINT (7.2462 47.14107) |
| 117 | 5611649752 | coffee | Pelluch GmbH Kaffeemaschienen | POINT (7.54718 47.54892) |
118 rows × 4 columns
Show selected features on map¶
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=9,
tiles='EsriWorldTopoMap')
# Map settings
folium.GeoJson(
gdf,
name='map',
popup=folium.GeoJsonPopup(fields=['name', 'shop'])
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Same but for dairy products¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query
sql = """SELECT
h.osm_id,
h.shop,
h.name,
ST_Transform(h.way, 4326) AS geom
FROM planet_osm_point h
WHERE h.shop = 'dairy';"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)
# Dispose the engine
engine.dispose()
# Print the GeoDataFrame
gdf
| osm_id | shop | name | geom | |
|---|---|---|---|---|
| 0 | 3086652026 | dairy | Käserei Sax | POINT (9.45643 47.2294) |
| 1 | 6331614596 | dairy | Nossa Caschareia | POINT (9.59642 46.59588) |
| 2 | 7796298536 | dairy | Laiterie | POINT (6.48589 46.86461) |
| 3 | 12090655227 | dairy | Laiterie | POINT (6.89998 46.52691) |
| 4 | 9628876022 | dairy | Laiterie Modèle | POINT (7.01441 46.25005) |
| ... | ... | ... | ... | ... |
| 76 | 980819667 | dairy | Käserei Büel | POINT (7.28702 46.74436) |
| 77 | 4307193372 | dairy | Laiterie - Fromagerie Pasquier | POINT (7.06689 46.76837) |
| 78 | 1512803412 | dairy | None | POINT (7.1554 47.36553) |
| 79 | 10298620447 | dairy | Dorfkäserei Thundorf | POINT (8.96682 47.54753) |
| 80 | 2990133771 | dairy | Chäs-Hütte | POINT (9.01833 47.5641) |
81 rows × 4 columns
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=9,
tiles='CartoDBDarkMatter')
# Map settings
folium.GeoJson(
gdf,
name='map',
popup=folium.GeoJsonPopup(fields=['name', 'shop'])
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Query: Select all supermarkets in a distance of 1000m around the central station in the city of Winterthur.¶
Note:
For each supermarket, the distance to the central station in meters is calculated and stored as new column 'distance_meters'.
In addition, a popup field was added to the map, allowing users to view detailed information about each selected feature when they click on it.
The WGS84 (World Geodetic System 1984) coordinates in ST_MakePoint(LON, LAT) were derived from: https://tools.retorte.ch/map.
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query
sql = """SELECT
p.osm_id,
p.shop,
p.name,
ST_Distance(
ST_Transform(p.way, 4326)::geography,
-- Central station coordinates
ST_SetSRID(ST_MakePoint(8.72397, 47.50031), 4326)::geography
) AS distance_meters,
ST_TRANSFORM(p.way, 4326) AS geom
FROM
planet_osm_point AS p
WHERE
p.shop = 'supermarket'
AND ST_DWithin(
ST_Transform(p.way, 4326)::geography,
-- Central station coordinates
ST_SetSRID(ST_MakePoint(8.72397, 47.50031), 4326)::geography,
1000
)
ORDER BY distance_meters;"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)
# Dispose the engine
engine.dispose()
# Print the GeoDataFrame
gdf
| osm_id | shop | name | distance_meters | geom | |
|---|---|---|---|---|---|
| 0 | 706203439 | supermarket | Coop | 159.883419 | POINT (8.72594 47.50085) |
| 1 | 4109460421 | supermarket | Asia Shop | 162.391281 | POINT (8.72208 47.50101) |
| 2 | 3831772784 | supermarket | Migros | 247.578208 | POINT (8.72115 47.49916) |
| 3 | 7380954145 | supermarket | Alnatura | 256.838011 | POINT (8.72074 47.49958) |
| 4 | 4095400190 | supermarket | ALDI | 274.275393 | POINT (8.72476 47.4979) |
| 5 | 4125136758 | supermarket | Tandoor Indischer Supermarkt | 290.212664 | POINT (8.72017 47.50073) |
| 6 | 4095400136 | supermarket | Denner | 316.354037 | POINT (8.72036 47.49886) |
| 7 | 709022324 | supermarket | Claro Weltladen | 441.129317 | POINT (8.72912 47.49842) |
| 8 | 4058248551 | supermarket | Migros | 600.117307 | POINT (8.73193 47.50012) |
| 9 | 3441033104 | supermarket | L'Ultimo Bacio | 680.202961 | POINT (8.73299 47.49999) |
Show selected features on map¶
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=16,
tiles='ESRIWorldImagery')
# Map settings
folium.GeoJson(
gdf,
name='map',
popup=folium.GeoJsonPopup(fields=['name', 'distance_meters'])
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Same but for coffee aroudn my home¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query
sql = """SELECT
p.osm_id,
p.shop,
p.name,
ST_Distance(
ST_Transform(p.way, 4326)::geography,
-- my coordinates
ST_SetSRID(ST_MakePoint(8.522117, 47.370249), 4326)::geography
) AS distance_meters,
ST_TRANSFORM(p.way, 4326) AS geom
FROM
planet_osm_point AS p
WHERE
p.shop = 'coffee'
AND ST_DWithin(
ST_Transform(p.way, 4326)::geography,
-- my coordinates
ST_SetSRID(ST_MakePoint(8.522117, 47.370249), 4326)::geography,
2000
)
ORDER BY distance_meters;"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine)
# Dispose the engine
engine.dispose()
# Print the GeoDataFrame
gdf
| osm_id | shop | name | distance_meters | geom | |
|---|---|---|---|---|---|
| 0 | 4370365690 | coffee | Kaffeepur | 277.208683 | POINT (8.51922 47.37178) |
| 1 | 4515947589 | coffee | Kaffeewerkstatt | 351.997864 | POINT (8.5175 47.3707) |
| 2 | 5297627699 | coffee | Rent a Barista | 797.665201 | POINT (8.53152 47.37352) |
| 3 | 4891858465 | coffee | Tchibo | 1017.406427 | POINT (8.535 47.36757) |
| 4 | 730380415 | coffee | ViCAFE | 1053.386237 | POINT (8.53551 47.3676) |
| 5 | 10967856734 | coffee | Beanbank | 1093.757918 | POINT (8.53654 47.36939) |
| 6 | 10795431343 | coffee | Bean Bank Coffee & CO | 1182.673128 | POINT (8.53304 47.37787) |
| 7 | 11350560294 | coffee | Nespresso Boutique | 1309.342540 | POINT (8.53812 47.37478) |
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=16,
tiles='ESRIWorldImagery')
# Map settings
folium.GeoJson(
gdf,
name='map',
popup=folium.GeoJsonPopup(fields=['name', 'distance_meters'])
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Query: Select all roads classified as 'motorway' and create a 5000m buffer around these roads.¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query (major roads)
sql = """-- Create buffer around major roads
SELECT
1 as group_id,
ST_TRANSFORM(ST_UNION(ST_Buffer(p.way::geometry, 5000)), 4326) AS geom
FROM public.planet_osm_roads AS p
WHERE
highway = 'motorway';"""
# Query the database and store the result in a GeoDataFrame
gdf = gpd.GeoDataFrame.from_postgis(sql, engine, geom_col='geom')
# Dispose the engine
engine.dispose()
Show selected features on map¶
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if gdf.crs is None:
gdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = gdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=9,
tiles='EsriWorldTopoMap')
# Map settings
folium.GeoJson(
gdf,
name='map'
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Same ercise but with other type of roads.¶
# Create SQLAlchemy Engine
engine = create_engine(db_connection_url)
# Define SQL query (major roads)
sql = """-- Create buffer around major roads
SELECT
1 as group_id,
ST_TRANSFORM(ST_UNION(ST_Buffer(p.way::geometry, 2000)), 4326) AS geom
FROM public.planet_osm_roads AS p
WHERE
highway = 'primary';"""
# Query the database and store the result in a GeoDataFrame
primarygdf = gpd.GeoDataFrame.from_postgis(sql, engine, geom_col='geom')
# Dispose the engine
engine.dispose()
# Ensure the GeoDataFrame has the correct projection (EPSG:4326)
if primarygdf.crs is None:
primarygdf.set_crs(epsg=4326, inplace=True)
else:
pass
# Calculate the mean longitude and latitude for the map center
centroids = primarygdf.geometry.centroid
lon = centroids.x.mean()
lat = centroids.y.mean()
# Initialize the map
m = folium.Map(location=[lat, lon],
zoom_start=9,
tiles='EsriWorldTopoMap')
folium.GeoJson(
primarygdf,
name='map',
fill_color='green',
fillOpacity=0.4,
).add_to(m)
folium.GeoJson(
gdf,
name='map',
fill_color='red',
fillOpacity=0.4,
).add_to(m)
folium.LayerControl().add_to(m)
# Plot map
m
Jupyter notebook --footer info-- (please always provide this at the end of each notebook)¶
import os
import platform
import socket
from platform import python_version
from datetime import datetime
print('-----------------------------------')
print(os.name.upper())
print(platform.system(), '|', platform.release())
print('Datetime:', datetime.now().strftime("%Y-%m-%d %H:%M:%S"))
print('Python Version:', python_version())
print('-----------------------------------')
----------------------------------- POSIX Linux | 6.5.0-1025-azure Datetime: 2024-10-06 13:34:56 Python Version: 3.12.1 -----------------------------------